Model Fitting and Diagnostics

Packages

Code
# Ensures the package "pacman" is installed
if (!require("pacman")) install.packages("pacman")
Loading required package: pacman
Code
### Install/Load packages
pacman::p_load(dplyr,
               tidyr,
               brms,
               bayesplot,
               tidybayes,
               ggplot2,
               here,
               MetBrewer,
               stringr)

color_scheme_set("red")

theme_set(
  theme(
    plot.title.position = 'plot',
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(size = 12),
    axis.text.y = element_text(size = 12),
    axis.title.x = element_text(size = 14),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", linewidth = 0.3)
  )
)

Data

Code
data = read.csv(here("data/cef_dataonly.csv"), header=T)

# Risk of bias data
rob = read.csv(here("data/ROB_cefepime.csv"), header=T)

filtered_data = 
  data |> 
  left_join(rob, by = "Study.ID") |> 
  filter(Comparator != "Control", #exclude sponsor studies
         N_cef != "NA") |>  # exclude studies with missing data
  rename("death_com" = deaths_com,
         "study" = Study.ID,
         "indication" = Intended.clinical.indication,
         "comparator_group" = Comparator.group,
         "population" = Patient.population,
         "overall_rob" = Overall.ROB) 

long_data = 
  filtered_data |> 
  pivot_longer(
   cols = c(death_cef, death_com, N_cef, N_com),
   names_to = c(".value", "group"),
   names_sep = "_"
) |> 
  # important for model fitting
  mutate(study = factor(study),
         treat = ifelse(group == "com", 0, 1) |> factor(levels = c("0", "1")),
         treat12 = ifelse(group == "com", -0.5, 0.5),
         
         indication = factor(indication,
                             levels = c("Febrile neutropenia",
                                        "Pneumonia",
                                        "Severe bacterial infections",
                                        "UTI",
                                        "Meningitis",      
                                        "Mixed")),
         comparator_group = factor(comparator_group,
                                   levels = c("Piperacillin-tazobactam",
                                              "Carbapenem",
                                              "Cefotaxime/Ceftriaxone",
                                              "Ceftazidime")),
         population = factor(population,
                             levels = c("Adults",
                                        "Pediatrics")),
         Cefepime_q12 = factor(Cefepime_q12,
                               levels = c("High",
                                          "Low")),
         
         overall_rob = factor(overall_rob,
                              levels = c("Low",
                                         "Some concerns",
                                         "High"))
         )

Overall Model

Model formula

Code
# Model 4 in DOI: 10.1002/sim.7588
mf = 
  bf(death | trials(N) ~ 1 + study + treat + (treat12 - 1 | study))

Priors

Code
get_prior(mf, family = binomial, data = long_data)
                prior     class                       coef group resp dpar
               (flat)         b                                           
               (flat)         b              studyAoun1997                
               (flat)         b              studyAtay2004                
               (flat)         b           studyAufiero1997                
               (flat)         b     studyAvilesMRobles2019                
               (flat)         b           studyBarckow1993                
               (flat)         b         studyBeaucaire1999                
               (flat)         b             studyBiron1998                
               (flat)         b             studyBohme1998                
               (flat)         b          studyBonfitto1999                
               (flat)         b               studyBow2006                
               (flat)         b           studyBradley2019                
               (flat)         b             studyBreen1997                
               (flat)         b         studyCannavino2015                
               (flat)         b      studyChandrasekar2000                
               (flat)         b             studyChang1998                
               (flat)         b            studyCherif2004                
               (flat)         b            studyChuang2002                
               (flat)         b           studyCordero2001                
               (flat)         b        studyCordonnier1997                
               (flat)         b           studyCornely2001                
               (flat)         b           studyCornely2002                
               (flat)         b         studyEdelstein1991                
               (flat)         b             studyErman2001                
               (flat)         b            studyFujita2016                
               (flat)         b            studyGentry1991                
               (flat)         b            studyGentry1992                
               (flat)         b           studyGhalaut2007                
               (flat)         b           studyGlauser1997                
               (flat)         b             studyGomez2001                
               (flat)         b             studyGomez2010                
               (flat)         b          studyGrossman1999                
               (flat)         b         studyHoepelman1993                
               (flat)         b          studyHolloway1996                
               (flat)         b             studyHuang2002                
               (flat)         b             studyHuang2005                
               (flat)         b              studyJehn1998                
               (flat)         b             studyJiang2003                
               (flat)         b            studyJindal2016                
               (flat)         b            studyKebudi2001                
               (flat)         b             studyKieft1994                
               (flat)         b            studyKutluk2004                
               (flat)         b              studyKwon2008                
               (flat)         b               studyLee2018                
               (flat)         b         studyLeophonte1993                
               (flat)         b               studyLin2001                
               (flat)         b               studyLiu2019                
               (flat)         b            studyMallet1997                
               (flat)         b           studyMcCabe1996a                
               (flat)         b           studyMcCabe1996b                
               (flat)         b           studyMustafa2001                
               (flat)         b           studyNakane2015a                
               (flat)         b           studyNakane2015b                
               (flat)         b            studyNaseem2011                
               (flat)         b            studyNewton1993                
               (flat)         b              studyOguz2006                
               (flat)         b                studyOi2020                
               (flat)         b             studyORyan1996                
               (flat)         b          studyPaladino2007                
               (flat)         b     studyPonceMdeMLeon1999                
               (flat)         b           studyPreheim1995                
               (flat)         b              studyQian2023                
               (flat)         b              studyRaad2003                
               (flat)         b           studyRamphal1997                
               (flat)         b      studySaezMLlorens1995                
               (flat)         b      studySaezMLlorens2001                
               (flat)         b            studySaito1992a                
               (flat)         b            studySaito1992b                
               (flat)         b              studySanz2002                
               (flat)         b         studySarashina2014                
               (flat)         b            studySchaad1998                
               (flat)         b           studySchrank1995                
               (flat)         b          studySchwartz1996                
               (flat)         b              studySeo2017a                
               (flat)         b              studySeo2017b                
               (flat)         b studyShahid2008DCPM4495001                
               (flat)         b           studySharifi1996                
               (flat)         b            studyTalaie2008                
               (flat)         b            studyTamura2002                
               (flat)         b             studyUygun2009                
               (flat)         b              studyWang1999                
               (flat)         b            studyWillis1998                
               (flat)         b          studyYakovlev2006                
               (flat)         b           studyZanetti2003                
               (flat)         b            studyZervos1998                
               (flat)         b                     treat1                
 student_t(3, 0, 2.5) Intercept                                           
 student_t(3, 0, 2.5)        sd                                           
 student_t(3, 0, 2.5)        sd                            study          
 student_t(3, 0, 2.5)        sd                    treat12 study          
 nlpar lb ub       source
                  default
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
                  default
        0         default
        0    (vectorized)
        0    (vectorized)
Code
informative = bayesmeta::TurnerEtAlPrior("all-cause mortality",
                                         "pharma",
                                         "pharma")


logmean = informative$parameters["tau", "mu"]
logsd = informative$parameters["tau", "sigma"]

# logmean
# -2.09

# logsd
# 0.705

priors_overall = 
  # study's baseline (treat = 0) log odds
   prior(normal(0, 1.5), class = Intercept) + 
  prior(normal(0, 1.5), class = b) + 
  # overall treatment effect, log odds ratio, "theta"
  prior(normal(0, 0.82), class = b, coef = "treat1") + 
  # between-study SD, "tau"
  prior(lognormal(-2.09, 0.705), class = sd)

Fit the model only sampling from prior distributions

Code
overall_mod_prior = 
  brm(
    sample_prior = "only",
    formula = mf,
    family = binomial,
    prior = priors_overall, 
    data = long_data,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "overall_mod_prior.Rds"),
    file_refit = "on_change"
)

Now check if the estimand distribution looks reasonable (prior predictive check):

Code
as_draws_df(overall_mod_prior) |> 
  ggplot() +
  aes(x = b_treat1) +
  stat_halfeye(.width = 0.95,
               point_interval = median_hdi) + 
  scale_x_continuous(
    breaks = c(0.1, 0.2, 0.5, 1, 2, 5, 10) |> log(),
    labels = function(x) exp(x)
    ) +
  coord_cartesian(x = c(0.1, 10) |> log())  +
  labs(x = "\nOdds Ratio (log scale)",
       y = " ",
       title = "Prior Predictive Check of the Treatment Effect Estimand") +
  theme(axis.text.y = element_blank())
Loading required package: rstan
Loading required package: StanHeaders

rstan version 2.35.0.9000 (Stan version 2.35.0)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
change `threads_per_chain` option:
rstan_options(threads_per_chain = 1)

Attaching package: 'rstan'
The following object is masked from 'package:tidyr':

    extract

Fit the full model (data + prior)

Code
overall_mod = 
  brm(
    formula = mf,
    family = binomial,
    prior = priors_overall, 
    data = long_data,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "overall_mod.Rds"),
    file_refit = "on_change"
)

Model Diagnostics

Rhat

Code
rhat(overall_mod$fit) |> 
  mcmc_rhat_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios of effective sample size to total sample size

Values are colored using different shades (lighter is better). The chosen thresholds are somewhat arbitrary, but can be useful guidelines in practice.

light: between 0.5 and 1 (high)

mid: between 0.1 and 0.5 (good)

dark: below 0.1 (low)

Code
neff_ratio(overall_mod) |> 
  mcmc_neff_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Does the Trace-Plot Exhibit Convergence?

Code
mcmc_trace(overall_mod, pars = c("b_treat1","sd_study__treat12"))

Does the Histogram Have Enough Information?

Code
mcmc_hist(overall_mod, pars = c("b_treat1","sd_study__treat12"))
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Do the Chains Exhibit a Strong Degree of Autocorrelation?

Code
mcmc_acf(overall_mod, pars = c("b_treat1","sd_study__treat12"),
         lags = 10)

Does the Posterior Distribution Make Substantive Sense?

Code
mcmc_dens(overall_mod,
          pars = c("b_treat1","sd_study__treat12"),
          transformations = list("b_treat1" = "exp"))

Meta-Regressions

Clinical Indication

Model formula

Code
# Formula with meta-regression
mf_indication = 
  bf(death | trials(N) ~ 1 + study + treat*indication + (treat12 - 1 | study))

get_prior(mf_indication, data = long_data, family = binomial)
                prior     class                                       coef
               (flat)         b                                           
               (flat)         b                       indicationMeningitis
               (flat)         b                            indicationMixed
               (flat)         b                        indicationPneumonia
               (flat)         b        indicationSeverebacterialinfections
               (flat)         b                              indicationUTI
               (flat)         b                              studyAoun1997
               (flat)         b                              studyAtay2004
               (flat)         b                           studyAufiero1997
               (flat)         b                     studyAvilesMRobles2019
               (flat)         b                           studyBarckow1993
               (flat)         b                         studyBeaucaire1999
               (flat)         b                             studyBiron1998
               (flat)         b                             studyBohme1998
               (flat)         b                          studyBonfitto1999
               (flat)         b                               studyBow2006
               (flat)         b                           studyBradley2019
               (flat)         b                             studyBreen1997
               (flat)         b                         studyCannavino2015
               (flat)         b                      studyChandrasekar2000
               (flat)         b                             studyChang1998
               (flat)         b                            studyCherif2004
               (flat)         b                            studyChuang2002
               (flat)         b                           studyCordero2001
               (flat)         b                        studyCordonnier1997
               (flat)         b                           studyCornely2001
               (flat)         b                           studyCornely2002
               (flat)         b                         studyEdelstein1991
               (flat)         b                             studyErman2001
               (flat)         b                            studyFujita2016
               (flat)         b                            studyGentry1991
               (flat)         b                            studyGentry1992
               (flat)         b                           studyGhalaut2007
               (flat)         b                           studyGlauser1997
               (flat)         b                             studyGomez2001
               (flat)         b                             studyGomez2010
               (flat)         b                          studyGrossman1999
               (flat)         b                         studyHoepelman1993
               (flat)         b                          studyHolloway1996
               (flat)         b                             studyHuang2002
               (flat)         b                             studyHuang2005
               (flat)         b                              studyJehn1998
               (flat)         b                             studyJiang2003
               (flat)         b                            studyJindal2016
               (flat)         b                            studyKebudi2001
               (flat)         b                             studyKieft1994
               (flat)         b                            studyKutluk2004
               (flat)         b                              studyKwon2008
               (flat)         b                               studyLee2018
               (flat)         b                         studyLeophonte1993
               (flat)         b                               studyLin2001
               (flat)         b                               studyLiu2019
               (flat)         b                            studyMallet1997
               (flat)         b                           studyMcCabe1996a
               (flat)         b                           studyMcCabe1996b
               (flat)         b                           studyMustafa2001
               (flat)         b                           studyNakane2015a
               (flat)         b                           studyNakane2015b
               (flat)         b                            studyNaseem2011
               (flat)         b                            studyNewton1993
               (flat)         b                              studyOguz2006
               (flat)         b                                studyOi2020
               (flat)         b                             studyORyan1996
               (flat)         b                          studyPaladino2007
               (flat)         b                     studyPonceMdeMLeon1999
               (flat)         b                           studyPreheim1995
               (flat)         b                              studyQian2023
               (flat)         b                              studyRaad2003
               (flat)         b                           studyRamphal1997
               (flat)         b                      studySaezMLlorens1995
               (flat)         b                      studySaezMLlorens2001
               (flat)         b                            studySaito1992a
               (flat)         b                            studySaito1992b
               (flat)         b                              studySanz2002
               (flat)         b                         studySarashina2014
               (flat)         b                            studySchaad1998
               (flat)         b                           studySchrank1995
               (flat)         b                          studySchwartz1996
               (flat)         b                              studySeo2017a
               (flat)         b                              studySeo2017b
               (flat)         b                 studyShahid2008DCPM4495001
               (flat)         b                           studySharifi1996
               (flat)         b                            studyTalaie2008
               (flat)         b                            studyTamura2002
               (flat)         b                             studyUygun2009
               (flat)         b                              studyWang1999
               (flat)         b                            studyWillis1998
               (flat)         b                          studyYakovlev2006
               (flat)         b                           studyZanetti2003
               (flat)         b                            studyZervos1998
               (flat)         b                                     treat1
               (flat)         b                treat1:indicationMeningitis
               (flat)         b                     treat1:indicationMixed
               (flat)         b                 treat1:indicationPneumonia
               (flat)         b treat1:indicationSeverebacterialinfections
               (flat)         b                       treat1:indicationUTI
 student_t(3, 0, 2.5) Intercept                                           
 student_t(3, 0, 2.5)        sd                                           
 student_t(3, 0, 2.5)        sd                                           
 student_t(3, 0, 2.5)        sd                                    treat12
 group resp dpar nlpar lb ub       source
                                  default
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                                  default
                        0         default
 study                  0    (vectorized)
 study                  0    (vectorized)

Priors

Code
priors_regression_indication = 
  # study's baseline (treat = 0) log odds, when subgroup == "Febrile neutropenia"
  prior(normal(0, 1.5), class = Intercept) +
  prior(normal(0, 1.5), class = b) + 
  # treatment effect in reference subgroup "Febrile neutropenia", log odds ratio
  prior(normal(0, 0.82), class = b, coef = "treat1") + 
  # change in log-odds when treat = 0 for each subgroup
  prior(normal(0, 1.2), class = b, coef = "indicationMeningitis") + 
  prior(normal(0, 1.2), class = b, coef = "indicationMixed") + 
  prior(normal(0, 1.2), class = b, coef = "indicationPneumonia") + 
  prior(normal(0, 1.2), class = b, coef = "indicationSeverebacterialinfections") + 
  prior(normal(0, 1.2), class = b, coef = "indicationUTI") + 
  # change in the treatment effect for each subgroup compared to "Febrile neutropenia"
  prior(normal(0, 0.3), class = b, coef = "treat1:indicationMeningitis") + 
  prior(normal(0, 0.3), class = b, coef = "treat1:indicationMixed") + 
  prior(normal(0, 0.3), class = b, coef = "treat1:indicationPneumonia") + 
  prior(normal(0, 0.3), class = b, coef = "treat1:indicationSeverebacterialinfections") + 
  prior(normal(0, 0.3), class = b, coef = "treat1:indicationUTI") + 
  # between-study SD, "tau"
  prior(lognormal(-2.09, 0.705), class = sd)

Fit the model only sampling from prior distributions

Code
indication_mod_prior = 
  brm(
    sample_prior = "only",
    formula = mf_indication,
    family = binomial,
    prior = priors_regression_indication, 
    data = long_data,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "indication_mod_prior.Rds"),
    file_refit = "on_change"
)

Now check if estimand distributions look reasonable (prior predictive check):

Code
as_draws_df(indication_mod_prior) |> 
  reframe("Febrile Neutropenia" = b_treat1,
          "Mixed" = b_treat1 + `b_treat1:indicationMixed`,
          "Pneumonia" = b_treat1 + `b_treat1:indicationPneumonia`,
          "Severe bacterial infections" = b_treat1 + `b_treat1:indicationSeverebacterialinfections`,
          "UTI" = b_treat1 + `b_treat1:indicationUTI`,
          "Meningitis" = b_treat1 + `b_treat1:indicationMeningitis`) |> 
  pivot_longer(everything()) |> 
  mutate(name = factor(name, levels = c("Febrile Neutropenia",                     
                                        "Pneumonia",
                                        "Severe bacterial infections",
                                        "UTI",  
                                        "Mixed",
                                        "Meningitis") |> rev()
                       )
         ) |> 
  
  ggplot() +
  aes(x = value, y = name, fill = name) +
  stat_halfeye(.width = 0.95,
               point_interval = median_hdi) +
  scale_fill_manual(values=met.brewer("Isfahan1", 8) |> rev()) +
  scale_x_continuous(
    breaks = c(0.1, 0.2, 0.5, 1, 2, 5, 10) |> log(),
    labels = function(x) exp(x)
    ) +
  coord_cartesian(x = c(0.1, 10) |> log())  +
  labs(x = "\nOdds Ratio (log scale)",
       y = " ") +
  theme(legend.position = "none",
        axis.text.y = element_text(hjust = 1))

Fit the full model (data + prior)

Code
indication_mod = 
  brm(
    formula = mf_indication,
    family = binomial,
    prior = priors_regression_indication, 
    data = long_data,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "indication_mod.Rds"),
    file_refit = "on_change"
)

Model Diagnostics

Rhat

Code
rhat(indication_mod$fit) |> 
  mcmc_rhat_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios of effective sample size to total sample size

Values are colored using different shades (lighter is better). The chosen thresholds are somewhat arbitrary, but can be useful guidelines in practice.

light: between 0.5 and 1 (high)

mid: between 0.1 and 0.5 (good)

dark: below 0.1 (low)

Code
neff_ratio(indication_mod) |> 
  mcmc_neff_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Does the Trace-Plot Exhibit Convergence?

Code
mcmc_trace(indication_mod, pars = c("b_treat1",
                                    "b_treat1:indicationMeningitis",
                                    "b_treat1:indicationMixed",
                                    "b_treat1:indicationPneumonia",
                                    "b_treat1:indicationSeverebacterialinfections",
                                    "b_treat1:indicationUTI",
                                    "sd_study__treat12"))

Does the Histogram Have Enough Information?

Code
mcmc_hist(indication_mod, pars = c("b_treat1",
                                    "b_treat1:indicationMeningitis",
                                    "b_treat1:indicationMixed",
                                    "b_treat1:indicationPneumonia",
                                    "b_treat1:indicationSeverebacterialinfections",
                                    "b_treat1:indicationUTI",
                                    "sd_study__treat12"))
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Do the Chains Exhibit a Strong Degree of Autocorrelation?

Code
mcmc_acf(indication_mod, pars = c("b_treat1",
                                  "b_treat1:indicationMeningitis",
                                  "b_treat1:indicationMixed",
                                  "b_treat1:indicationPneumonia",
                                  "b_treat1:indicationSeverebacterialinfections",
                                  "b_treat1:indicationUTI",
                                  "sd_study__treat12"),
         lags = 10)

Does the Posterior Distribution Make Substantive Sense?

Code
mcmc_dens(indication_mod,
          pars = c("b_treat1",
                   "sd_study__treat12"),
          transformations = list("b_treat1" = "exp"))

Comparator Group

Model formula

Code
# Formula with meta-regression
mf_comparator = 
  bf(death | trials(N) ~ 1 + study + treat*comparator_group + (treat12 - 1 | study))

get_prior(mf_comparator, data = long_data, family = binomial)
Warning: Rows containing NAs were excluded from the model.
                prior     class                                          coef
               (flat)         b                                              
               (flat)         b                    comparator_groupCarbapenem
               (flat)         b        comparator_groupCefotaximeDCeftriaxone
               (flat)         b                   comparator_groupCeftazidime
               (flat)         b                                 studyAoun1997
               (flat)         b                                 studyAtay2004
               (flat)         b                              studyAufiero1997
               (flat)         b                              studyBarckow1993
               (flat)         b                            studyBeaucaire1999
               (flat)         b                                studyBiron1998
               (flat)         b                                studyBohme1998
               (flat)         b                             studyBonfitto1999
               (flat)         b                                  studyBow2006
               (flat)         b                            studyCannavino2015
               (flat)         b                         studyChandrasekar2000
               (flat)         b                                studyChang1998
               (flat)         b                               studyCherif2004
               (flat)         b                               studyChuang2002
               (flat)         b                              studyCordero2001
               (flat)         b                           studyCordonnier1997
               (flat)         b                              studyCornely2001
               (flat)         b                              studyCornely2002
               (flat)         b                            studyEdelstein1991
               (flat)         b                                studyErman2001
               (flat)         b                               studyFujita2016
               (flat)         b                               studyGentry1991
               (flat)         b                               studyGentry1992
               (flat)         b                              studyGhalaut2007
               (flat)         b                              studyGlauser1997
               (flat)         b                                studyGomez2001
               (flat)         b                                studyGomez2010
               (flat)         b                             studyGrossman1999
               (flat)         b                            studyHoepelman1993
               (flat)         b                             studyHolloway1996
               (flat)         b                                studyHuang2002
               (flat)         b                                 studyJehn1998
               (flat)         b                                studyJiang2003
               (flat)         b                               studyJindal2016
               (flat)         b                               studyKebudi2001
               (flat)         b                                studyKieft1994
               (flat)         b                               studyKutluk2004
               (flat)         b                                 studyKwon2008
               (flat)         b                                  studyLee2018
               (flat)         b                            studyLeophonte1993
               (flat)         b                                  studyLin2001
               (flat)         b                               studyMallet1997
               (flat)         b                              studyMcCabe1996a
               (flat)         b                              studyMcCabe1996b
               (flat)         b                              studyMustafa2001
               (flat)         b                              studyNakane2015b
               (flat)         b                               studyNewton1993
               (flat)         b                                 studyOguz2006
               (flat)         b                                   studyOi2020
               (flat)         b                                studyORyan1996
               (flat)         b                             studyPaladino2007
               (flat)         b                        studyPonceMdeMLeon1999
               (flat)         b                              studyPreheim1995
               (flat)         b                                 studyQian2023
               (flat)         b                                 studyRaad2003
               (flat)         b                              studyRamphal1997
               (flat)         b                         studySaezMLlorens1995
               (flat)         b                         studySaezMLlorens2001
               (flat)         b                               studySaito1992a
               (flat)         b                               studySaito1992b
               (flat)         b                                 studySanz2002
               (flat)         b                               studySchaad1998
               (flat)         b                              studySchrank1995
               (flat)         b                             studySchwartz1996
               (flat)         b                                 studySeo2017a
               (flat)         b                                 studySeo2017b
               (flat)         b                    studyShahid2008DCPM4495001
               (flat)         b                              studySharifi1996
               (flat)         b                               studyTalaie2008
               (flat)         b                               studyTamura2002
               (flat)         b                                studyUygun2009
               (flat)         b                                 studyWang1999
               (flat)         b                               studyWillis1998
               (flat)         b                             studyYakovlev2006
               (flat)         b                              studyZanetti2003
               (flat)         b                               studyZervos1998
               (flat)         b                                        treat1
               (flat)         b             treat1:comparator_groupCarbapenem
               (flat)         b treat1:comparator_groupCefotaximeDCeftriaxone
               (flat)         b            treat1:comparator_groupCeftazidime
 student_t(3, 0, 2.5) Intercept                                              
 student_t(3, 0, 2.5)        sd                                              
 student_t(3, 0, 2.5)        sd                                              
 student_t(3, 0, 2.5)        sd                                       treat12
 group resp dpar nlpar lb ub       source
                                  default
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                             (vectorized)
                                  default
                        0         default
 study                  0    (vectorized)
 study                  0    (vectorized)

Priors

Code
priors_regression_comparator = 
  # study's baseline (treat = 0) log odds when subgroup == "Piperacillin-tazobactam"
  prior(normal(0, 1.5), class = Intercept) +
  prior(normal(0, 1.5), class = b) + 
  # treatment effect in reference subgroup "Piperacillin-tazobactam", log odds ratio
  prior(normal(0, 0.82), class = b, coef = "treat1") + 
  # change in log-odds when treat = 0 for each subgroup
  prior(normal(0, 1.2), class = b, coef = "comparator_groupCarbapenem") + 
  prior(normal(0, 1.2), class = b, coef = "comparator_groupCefotaximeDCeftriaxone") + 
  prior(normal(0, 1.2), class = b, coef = "comparator_groupCeftazidime") + 
  # change in the treatment effect for each subgroup compared to "Piperacillin-tazobactam"
  prior(normal(0, 0.3), class = b, coef = "treat1:comparator_groupCarbapenem") + 
  prior(normal(0, 0.3), class = b, coef = "treat1:comparator_groupCefotaximeDCeftriaxone") + 
  prior(normal(0, 0.3), class = b, coef = "treat1:comparator_groupCeftazidime") + 
  # between-study SD, "tau"
  prior(lognormal(-2.09, 0.705), class = sd)

Fit the model only sampling from prior distributions

Code
comparator_mod_prior = 
  brm(
    sample_prior = "only",
    formula = mf_comparator,
    family = binomial,
    prior = priors_regression_comparator, 
    data = subset(long_data, !is.na(comparator_group)), # remove NA
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "comparator_mod_prior.Rds"),
    file_refit = "on_change"
)

Now check if estimand distributions look reasonable (prior predictive check):

Code
as_draws_df(comparator_mod_prior) |> 
  reframe("Piperacillin-tazobactam" = b_treat1,
          "Carbapenem" = b_treat1 + `b_treat1:comparator_groupCarbapenem`,
          "Cefotaxime/Ceftriaxone" = b_treat1 + `b_treat1:comparator_groupCefotaximeDCeftriaxone`,
          "Ceftazidime" = b_treat1 + `b_treat1:comparator_groupCeftazidime`) |> 
  pivot_longer(everything()) |> 
  mutate(name = factor(name, c("Piperacillin-tazobactam",
                               "Carbapenem",
                               "Cefotaxime/Ceftriaxone",
                               "Ceftazidime") |> rev()
                       )
         ) |> 
  
  ggplot() +
  aes(x = value, y = name, fill = name) +
  stat_halfeye(.width = 0.95,
               point_interval = median_hdi) +
  scale_fill_manual(values=met.brewer("Isfahan1", 4) |> rev()) +
  scale_x_continuous(
    breaks = c(0.1, 0.2, 0.5, 1, 2, 5, 10) |> log(),
    labels = function(x) exp(x)
    ) +
  coord_cartesian(x = c(0.1, 10) |> log())  +
  labs(x = "\nOdds Ratio (log scale)",
       y = " ") +
  theme(legend.position = "none",
        axis.text.y = element_text(hjust = 1))

Fit the full model (data + prior)

Code
comparator_mod = 
  brm(
    formula = mf_comparator,
    family = binomial,
    prior = priors_regression_comparator, 
    data = subset(long_data, !is.na(comparator_group)), # remove NA
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "comparator_mod.Rds"),
    file_refit = "on_change"
)

Model Diagnostics

Rhat

Code
rhat(comparator_mod$fit) |> 
  mcmc_rhat_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios of effective sample size to total sample size

Values are colored using different shades (lighter is better). The chosen thresholds are somewhat arbitrary, but can be useful guidelines in practice.

light: between 0.5 and 1 (high)

mid: between 0.1 and 0.5 (good)

dark: below 0.1 (low)

Code
neff_ratio(comparator_mod) |> 
  mcmc_neff_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Does the Trace-Plot Exhibit Convergence?

Code
mcmc_trace(comparator_mod, pars = c("b_treat1",
                                    "b_treat1:comparator_groupCarbapenem",
                                    "b_treat1:comparator_groupCefotaximeDCeftriaxone",
                                    "b_treat1:comparator_groupCeftazidime",
                                    "sd_study__treat12"))

Does the Histogram Have Enough Information?

Code
mcmc_hist(comparator_mod, pars = c("b_treat1",
                                    "b_treat1:comparator_groupCarbapenem",
                                    "b_treat1:comparator_groupCefotaximeDCeftriaxone",
                                    "b_treat1:comparator_groupCeftazidime",
                                    "sd_study__treat12"))
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Do the Chains Exhibit a Strong Degree of Autocorrelation?

Code
mcmc_acf(comparator_mod, pars = c("b_treat1",
                                    "b_treat1:comparator_groupCarbapenem",
                                    "b_treat1:comparator_groupCefotaximeDCeftriaxone",
                                    "b_treat1:comparator_groupCeftazidime",
                                    "sd_study__treat12"),
         lags = 10)

Does the Posterior Distribution Make Substantive Sense?

Code
mcmc_dens(comparator_mod,
          pars = c("b_treat1",
                   "sd_study__treat12"),
          transformations = list("b_treat1" = "exp"))

Population

Model formula

Code
# Formula with meta-regression
mf_population = 
  bf(death | trials(N) ~ 1 + study + treat*population + (treat12 - 1 | study))

get_prior(mf_population, data = long_data, family = binomial)
                prior     class                        coef group resp dpar
               (flat)         b                                            
               (flat)         b        populationPediatrics                
               (flat)         b               studyAoun1997                
               (flat)         b               studyAtay2004                
               (flat)         b            studyAufiero1997                
               (flat)         b      studyAvilesMRobles2019                
               (flat)         b            studyBarckow1993                
               (flat)         b          studyBeaucaire1999                
               (flat)         b              studyBiron1998                
               (flat)         b              studyBohme1998                
               (flat)         b           studyBonfitto1999                
               (flat)         b                studyBow2006                
               (flat)         b            studyBradley2019                
               (flat)         b              studyBreen1997                
               (flat)         b          studyCannavino2015                
               (flat)         b       studyChandrasekar2000                
               (flat)         b              studyChang1998                
               (flat)         b             studyCherif2004                
               (flat)         b             studyChuang2002                
               (flat)         b            studyCordero2001                
               (flat)         b         studyCordonnier1997                
               (flat)         b            studyCornely2001                
               (flat)         b            studyCornely2002                
               (flat)         b          studyEdelstein1991                
               (flat)         b              studyErman2001                
               (flat)         b             studyFujita2016                
               (flat)         b             studyGentry1991                
               (flat)         b             studyGentry1992                
               (flat)         b            studyGhalaut2007                
               (flat)         b            studyGlauser1997                
               (flat)         b              studyGomez2001                
               (flat)         b              studyGomez2010                
               (flat)         b           studyGrossman1999                
               (flat)         b          studyHoepelman1993                
               (flat)         b           studyHolloway1996                
               (flat)         b              studyHuang2002                
               (flat)         b              studyHuang2005                
               (flat)         b               studyJehn1998                
               (flat)         b              studyJiang2003                
               (flat)         b             studyJindal2016                
               (flat)         b             studyKebudi2001                
               (flat)         b              studyKieft1994                
               (flat)         b             studyKutluk2004                
               (flat)         b               studyKwon2008                
               (flat)         b                studyLee2018                
               (flat)         b          studyLeophonte1993                
               (flat)         b                studyLin2001                
               (flat)         b                studyLiu2019                
               (flat)         b             studyMallet1997                
               (flat)         b            studyMcCabe1996a                
               (flat)         b            studyMcCabe1996b                
               (flat)         b            studyMustafa2001                
               (flat)         b            studyNakane2015a                
               (flat)         b            studyNakane2015b                
               (flat)         b             studyNaseem2011                
               (flat)         b             studyNewton1993                
               (flat)         b               studyOguz2006                
               (flat)         b                 studyOi2020                
               (flat)         b              studyORyan1996                
               (flat)         b           studyPaladino2007                
               (flat)         b      studyPonceMdeMLeon1999                
               (flat)         b            studyPreheim1995                
               (flat)         b               studyQian2023                
               (flat)         b               studyRaad2003                
               (flat)         b            studyRamphal1997                
               (flat)         b       studySaezMLlorens1995                
               (flat)         b       studySaezMLlorens2001                
               (flat)         b             studySaito1992a                
               (flat)         b             studySaito1992b                
               (flat)         b               studySanz2002                
               (flat)         b          studySarashina2014                
               (flat)         b             studySchaad1998                
               (flat)         b            studySchrank1995                
               (flat)         b           studySchwartz1996                
               (flat)         b               studySeo2017a                
               (flat)         b               studySeo2017b                
               (flat)         b  studyShahid2008DCPM4495001                
               (flat)         b            studySharifi1996                
               (flat)         b             studyTalaie2008                
               (flat)         b             studyTamura2002                
               (flat)         b              studyUygun2009                
               (flat)         b               studyWang1999                
               (flat)         b             studyWillis1998                
               (flat)         b           studyYakovlev2006                
               (flat)         b            studyZanetti2003                
               (flat)         b             studyZervos1998                
               (flat)         b                      treat1                
               (flat)         b treat1:populationPediatrics                
 student_t(3, 0, 2.5) Intercept                                            
 student_t(3, 0, 2.5)        sd                                            
 student_t(3, 0, 2.5)        sd                             study          
 student_t(3, 0, 2.5)        sd                     treat12 study          
 nlpar lb ub       source
                  default
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
                  default
        0         default
        0    (vectorized)
        0    (vectorized)

Priors

Code
priors_regression_population = 
  # study's baseline (treat = 0) log odds when subgroup == "Adults"
  prior(normal(0, 1.5), class = Intercept) +
  prior(normal(0, 1.5), class = b) + 
  # treatment effect in reference subgroup "Adults", log odds ratio
  prior(normal(0, 0.82), class = b, coef = "treat1") + 
  # change in log-odds when treat = 0 for each subgroup
  prior(normal(0, 1.2), class = b, coef = "populationPediatrics") + 
  # change in the treatment effect for the "Pediatrics" subgroup
  prior(normal(0, 0.3), class = b, coef = "treat1:populationPediatrics") + 
  # between-study SD, "tau"
  prior(lognormal(-2.09, 0.705), class = sd)

Fit the model only sampling from prior distributions

Code
population_mod_prior = 
  brm(
    sample_prior = "only",
    formula = mf_population,
    family = binomial,
    prior = priors_regression_population, 
    data = long_data,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "population_mod_prior.Rds"),
    file_refit = "on_change"
)

Now check if estimand distributions look reasonable (prior predictive check):

Code
as_draws_df(population_mod_prior) |> 
  reframe("Adults" = b_treat1,
          "Pediatrics" = b_treat1 + `b_treat1:populationPediatrics`) |> 
  pivot_longer(everything()) |> 
  mutate(name = factor(name, c("Adults",
                               "Pediatrics") |> rev()
                       )
         ) |> 
  
  ggplot() +
  aes(x = value, y = name, fill = name) +
  stat_halfeye(.width = 0.95,
               point_interval = median_hdi) +
  scale_fill_manual(values=met.brewer("Isfahan1", 2) |> rev()) +
  scale_x_continuous(
    breaks = c(0.1, 0.2, 0.5, 1, 2, 5, 10) |> log(),
    labels = function(x) exp(x)
    ) +
  coord_cartesian(x = c(0.1, 10) |> log())  +
  labs(x = "\nOdds Ratio (log scale)",
       y = " ") +
  theme(legend.position = "none",
        axis.text.y = element_text(hjust = 1))

Fit the full model (data + prior)

Code
population_mod = 
  brm(
    formula = mf_population,
    family = binomial,
    prior = priors_regression_population, 
    data = long_data,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "population_mod.Rds"),
    file_refit = "on_change"
)

Model Diagnostics

Rhat

Code
rhat(population_mod$fit) |> 
  mcmc_rhat_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios of effective sample size to total sample size

Values are colored using different shades (lighter is better). The chosen thresholds are somewhat arbitrary, but can be useful guidelines in practice.

light: between 0.5 and 1 (high)

mid: between 0.1 and 0.5 (good)

dark: below 0.1 (low)

Code
neff_ratio(population_mod) |> 
  mcmc_neff_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Does the Trace-Plot Exhibit Convergence?

Code
mcmc_trace(population_mod, pars = c("b_treat1",
                                    "b_treat1:populationPediatrics",
                                    "sd_study__treat12"))

Does the Histogram Have Enough Information?

Code
mcmc_hist(population_mod, pars = c("b_treat1",
                                    "b_treat1:populationPediatrics",
                                    "sd_study__treat12"))
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Do the Chains Exhibit a Strong Degree of Autocorrelation?

Code
mcmc_acf(population_mod, pars = c("b_treat1",
                                    "b_treat1:populationPediatrics",
                                    "sd_study__treat12"),
         lags = 10)

Does the Posterior Distribution Make Substantive Sense?

Code
mcmc_dens(population_mod,
          pars = c("b_treat1",
                   "b_treat1:populationPediatrics",
                   "sd_study__treat12"),
          transformations = list("b_treat1" = "exp"))

Cefepime Dose

Model formula

Code
# Formula with meta-regression
mf_cefdose = 
  bf(death | trials(N) ~ 1 + study + treat*Cefepime_q12 + (treat12 - 1 | study))

get_prior(mf_cefdose, data = long_data, family = binomial)
Warning: Rows containing NAs were excluded from the model.
                prior     class                   coef group resp dpar nlpar lb
               (flat)         b                                                
               (flat)         b        Cefepime_q12Low                         
               (flat)         b       studyAufiero1997                         
               (flat)         b       studyBarckow1993                         
               (flat)         b     studyBeaucaire1999                         
               (flat)         b         studyBiron1998                         
               (flat)         b         studyBohme1998                         
               (flat)         b      studyBonfitto1999                         
               (flat)         b           studyBow2006                         
               (flat)         b       studyBradley2019                         
               (flat)         b         studyBreen1997                         
               (flat)         b  studyChandrasekar2000                         
               (flat)         b         studyChang1998                         
               (flat)         b        studyCherif2004                         
               (flat)         b       studyCordero2001                         
               (flat)         b    studyCordonnier1997                         
               (flat)         b       studyCornely2001                         
               (flat)         b       studyCornely2002                         
               (flat)         b     studyEdelstein1991                         
               (flat)         b         studyErman2001                         
               (flat)         b        studyFujita2016                         
               (flat)         b        studyGentry1991                         
               (flat)         b       studyGlauser1997                         
               (flat)         b         studyGomez2001                         
               (flat)         b         studyGomez2010                         
               (flat)         b      studyGrossman1999                         
               (flat)         b     studyHoepelman1993                         
               (flat)         b      studyHolloway1996                         
               (flat)         b         studyHuang2002                         
               (flat)         b          studyJehn1998                         
               (flat)         b         studyJiang2003                         
               (flat)         b        studyJindal2016                         
               (flat)         b         studyKieft1994                         
               (flat)         b          studyKwon2008                         
               (flat)         b     studyLeophonte1993                         
               (flat)         b           studyLin2001                         
               (flat)         b           studyLiu2019                         
               (flat)         b       studyMcCabe1996a                         
               (flat)         b       studyMcCabe1996b                         
               (flat)         b       studyNakane2015a                         
               (flat)         b       studyNakane2015b                         
               (flat)         b        studyNaseem2011                         
               (flat)         b        studyNewton1993                         
               (flat)         b            studyOi2020                         
               (flat)         b      studyPaladino2007                         
               (flat)         b studyPonceMdeMLeon1999                         
               (flat)         b       studyPreheim1995                         
               (flat)         b          studyQian2023                         
               (flat)         b          studyRaad2003                         
               (flat)         b       studyRamphal1997                         
               (flat)         b        studySaito1992a                         
               (flat)         b        studySaito1992b                         
               (flat)         b          studySanz2002                         
               (flat)         b       studySchrank1995                         
               (flat)         b      studySchwartz1996                         
               (flat)         b          studySeo2017a                         
               (flat)         b          studySeo2017b                         
               (flat)         b       studySharifi1996                         
               (flat)         b        studyTalaie2008                         
               (flat)         b        studyTamura2002                         
               (flat)         b          studyWang1999                         
               (flat)         b        studyWillis1998                         
               (flat)         b      studyYakovlev2006                         
               (flat)         b       studyZanetti2003                         
               (flat)         b        studyZervos1998                         
               (flat)         b                 treat1                         
               (flat)         b treat1:Cefepime_q12Low                         
 student_t(3, 0, 2.5) Intercept                                                
 student_t(3, 0, 2.5)        sd                                               0
 student_t(3, 0, 2.5)        sd                        study                  0
 student_t(3, 0, 2.5)        sd                treat12 study                  0
 ub       source
         default
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
         default
         default
    (vectorized)
    (vectorized)

Priors

Code
priors_regression_cefdose = 
  # study's baseline (treat = 0) log odds when subgroup == "High"
  prior(normal(0, 1.5), class = Intercept) +
  prior(normal(0, 1.5), class = b) + 
  # treatment effect in reference subgroup "High", log odds ratio
  prior(normal(0, 0.82), class = b, coef = "treat1") + 
  # change in log-odds when treat = 0 for each subgroup
  prior(normal(0, 1.2), class = b, coef = "Cefepime_q12Low") + 
  # change in the treatment effect for the "Low" subgroup
  prior(normal(0, 0.3), class = b, coef = "treat1:Cefepime_q12Low") + 
  # between-study SD, "tau"
  prior(lognormal(-2.09, 0.705), class = sd)

Fit the model only sampling from prior distributions

Code
cefdose_mod_prior = 
  brm(
    sample_prior = "only",
    formula = mf_cefdose,
    family = binomial,
    prior = priors_regression_cefdose, 
    data = subset(long_data, !is.na(Cefepime_q12)), # remove NA,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "cefdose_mod_prior.Rds"),
    file_refit = "on_change"
)

Now check if estimand distributions look reasonable (prior predictive check):

Code
as_draws_df(cefdose_mod_prior) |> 
  reframe("High Dose" = b_treat1,
          "Low Dose" = b_treat1 + `b_treat1:Cefepime_q12Low`) |> 
  pivot_longer(everything()) |> 
  mutate(name = factor(name, c("High Dose",
                               "Low Dose") |> rev()
                       )
         ) |> 
  
  ggplot() +
  aes(x = value, y = name, fill = name) +
  stat_halfeye(.width = 0.95,
               point_interval = median_hdi) +
  scale_fill_manual(values=met.brewer("Isfahan1", 2) |> rev()) +
  scale_x_continuous(
    breaks = c(0.1, 0.2, 0.5, 1, 2, 5, 10) |> log(),
    labels = function(x) exp(x)
    ) +
  coord_cartesian(x = c(0.1, 10) |> log())  +
  labs(x = "\nOdds Ratio (log scale)",
       y = " ") +
  theme(legend.position = "none",
        axis.text.y = element_text(hjust = 1))

Fit the full model (data + prior)

Code
cefdose_mod = 
  brm(
    formula = mf_cefdose,
    family = binomial,
    prior = priors_regression_cefdose, 
    data = subset(long_data, !is.na(Cefepime_q12)), # remove NA,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "cefdose_mod.Rds"),
    file_refit = "on_change"
)

Model Diagnostics

Rhat

Code
rhat(cefdose_mod$fit) |> 
  mcmc_rhat_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios of effective sample size to total sample size

Values are colored using different shades (lighter is better). The chosen thresholds are somewhat arbitrary, but can be useful guidelines in practice.

light: between 0.5 and 1 (high)

mid: between 0.1 and 0.5 (good)

dark: below 0.1 (low)

Code
neff_ratio(cefdose_mod) |> 
  mcmc_neff_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

There is a parameter with low Neff/N. Check if most important parameters have high Neff/N:

Code
neff_ratio(cefdose_mod, pars = c("b_treat1",
                                    "b_treat1:Cefepime_q12Low",
                                    "sd_study__treat12")) |> 
  mcmc_neff() +
  theme(legend.position = "none" )

Does the Trace-Plot Exhibit Convergence?

Code
mcmc_trace(cefdose_mod, pars = c("b_treat1",
                                    "b_treat1:Cefepime_q12Low",
                                    "sd_study__treat12"))

Does the Histogram Have Enough Information?

Code
mcmc_hist(cefdose_mod, pars = c("b_treat1",
                                "b_treat1:Cefepime_q12Low",
                                "sd_study__treat12"))
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Do the Chains Exhibit a Strong Degree of Autocorrelation?

Code
mcmc_acf(cefdose_mod, pars = c("b_treat1",
                                "b_treat1:Cefepime_q12Low",
                                "sd_study__treat12"),
         lags = 10)

Does the Posterior Distribution Make Substantive Sense?

Code
mcmc_dens(cefdose_mod, pars = c("b_treat1",
                                "b_treat1:Cefepime_q12Low",
                                "sd_study__treat12"),
          transformations = list("b_treat1" = "exp"))

Risk of Bias

Model formula

Code
# Formula with meta-regression
mf_rob = 
  bf(death | trials(N) ~ 1 + study + treat*overall_rob + (treat12 - 1 | study))

get_prior(mf_rob, data = long_data, family = binomial)
Warning: Rows containing NAs were excluded from the model.
                prior     class                           coef group resp dpar
               (flat)         b                                               
               (flat)         b                overall_robHigh                
               (flat)         b        overall_robSomeconcerns                
               (flat)         b               studyAufiero1997                
               (flat)         b         studyAvilesMRobles2019                
               (flat)         b               studyBarckow1993                
               (flat)         b             studyBeaucaire1999                
               (flat)         b                 studyBiron1998                
               (flat)         b                 studyBohme1998                
               (flat)         b              studyBonfitto1999                
               (flat)         b                   studyBow2006                
               (flat)         b               studyBradley2019                
               (flat)         b                 studyBreen1997                
               (flat)         b             studyCannavino2015                
               (flat)         b          studyChandrasekar2000                
               (flat)         b                 studyChang1998                
               (flat)         b                studyCherif2004                
               (flat)         b                studyChuang2002                
               (flat)         b               studyCordero2001                
               (flat)         b            studyCordonnier1997                
               (flat)         b               studyCornely2002                
               (flat)         b             studyEdelstein1991                
               (flat)         b                 studyErman2001                
               (flat)         b                studyFujita2016                
               (flat)         b                studyGentry1991                
               (flat)         b               studyGhalaut2007                
               (flat)         b                 studyGomez2010                
               (flat)         b              studyGrossman1999                
               (flat)         b             studyHoepelman1993                
               (flat)         b              studyHolloway1996                
               (flat)         b                 studyHuang2002                
               (flat)         b                studyJindal2016                
               (flat)         b                studyKebudi2001                
               (flat)         b                 studyKieft1994                
               (flat)         b                studyKutluk2004                
               (flat)         b                  studyKwon2008                
               (flat)         b             studyLeophonte1993                
               (flat)         b                   studyLin2001                
               (flat)         b                   studyLiu2019                
               (flat)         b               studyMustafa2001                
               (flat)         b                studyNewton1993                
               (flat)         b                  studyOguz2006                
               (flat)         b                    studyOi2020                
               (flat)         b              studyPaladino2007                
               (flat)         b         studyPonceMdeMLeon1999                
               (flat)         b               studyPreheim1995                
               (flat)         b                  studyQian2023                
               (flat)         b                  studyRaad2003                
               (flat)         b          studySaezMLlorens1995                
               (flat)         b          studySaezMLlorens2001                
               (flat)         b                  studySanz2002                
               (flat)         b             studySarashina2014                
               (flat)         b                studySchaad1998                
               (flat)         b               studySchrank1995                
               (flat)         b              studySchwartz1996                
               (flat)         b               studySharifi1996                
               (flat)         b                studyTalaie2008                
               (flat)         b                studyTamura2002                
               (flat)         b                 studyUygun2009                
               (flat)         b                  studyWang1999                
               (flat)         b              studyYakovlev2006                
               (flat)         b               studyZanetti2003                
               (flat)         b                studyZervos1998                
               (flat)         b                         treat1                
               (flat)         b         treat1:overall_robHigh                
               (flat)         b treat1:overall_robSomeconcerns                
 student_t(3, 0, 2.5) Intercept                                               
 student_t(3, 0, 2.5)        sd                                               
 student_t(3, 0, 2.5)        sd                                study          
 student_t(3, 0, 2.5)        sd                        treat12 study          
 nlpar lb ub       source
                  default
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
                  default
        0         default
        0    (vectorized)
        0    (vectorized)

Priors

Code
priors_regression_rob = 
  # study's baseline (treat = 0) log odds when subgroup == "Low"
  prior(normal(0, 1.5), class = Intercept) +
  prior(normal(0, 1.5), class = b) + 
  # treatment effect in reference subgroup "Low", log odds ratio
  prior(normal(0, 0.82), class = b, coef = "treat1") + 
  # change in log-odds when treat = 0 for each subgroup
  prior(normal(0, 1.2), class = b, coef = "overall_robHigh") + 
  prior(normal(0, 1.2), class = b, coef = "overall_robSomeconcerns") + 
  # change in the treatment effect for each subgroup
  prior(normal(0, 0.3), class = b, coef = "treat1:overall_robHigh") + 
  prior(normal(0, 0.3), class = b, coef = "treat1:overall_robSomeconcerns") + 
  # between-study SD, "tau"
  prior(lognormal(-2.09, 0.705), class = sd)

Fit the model only sampling from prior distributions

Code
rob_mod_prior = 
  brm(
    sample_prior = "only",
    formula = mf_rob,
    family = binomial,
    prior = priors_regression_rob, 
    data = subset(long_data, !is.na(overall_rob)), # remove NA,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "rob_mod_prior.Rds"),
    file_refit = "on_change"
  )

Now check if estimand distributions look reasonable (prior predictive check):

Code
as_draws_df(rob_mod_prior) |> 
  reframe("Low" = b_treat1,
          "Some concerns" = b_treat1 + `b_treat1:overall_robSomeconcerns`,
          "High" = b_treat1 + `b_treat1:overall_robHigh`) |> 
  pivot_longer(everything()) |> 
  mutate(name = factor(name, c("Low",
                               "Some concerns",
                               "High") |> rev()
  )
  ) |> 
  
  ggplot() +
  aes(x = value, y = name, fill = name) +
  stat_halfeye(.width = 0.95,
               point_interval = median_hdi) +
  scale_fill_manual(values=met.brewer("Isfahan1", 3) |> rev()) +
  scale_x_continuous(
    breaks = c(0.1, 0.2, 0.5, 1, 2, 5, 10) |> log(),
    labels = function(x) exp(x)
  ) +
  coord_cartesian(x = c(0.1, 10) |> log())  +
  labs(x = "\nOdds Ratio (log scale)",
       y = " ") +
  theme(legend.position = "none",
        axis.text.y = element_text(hjust = 1))

Fit the full model (data + prior)

Code
rob_mod = 
  brm(
    formula = mf_rob,
    family = binomial,
    prior = priors_regression_rob, 
    data = subset(long_data, !is.na(overall_rob)), # remove NA,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "rob_mod.Rds"),
    file_refit = "on_change"
  )

Model Diagnostics

Rhat

Code
rhat(rob_mod$fit) |> 
  mcmc_rhat_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios of effective sample size to total sample size

Values are colored using different shades (lighter is better). The chosen thresholds are somewhat arbitrary, but can be useful guidelines in practice.

light: between 0.5 and 1 (high)

mid: between 0.1 and 0.5 (good)

dark: below 0.1 (low)

Code
neff_ratio(rob_mod) |> 
  mcmc_neff_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

There is a parameter with low Neff/N. Check if most important parameters have high Neff/N:

Code
neff_ratio(rob_mod, pars = c("b_treat1",
                             "b_treat1:overall_robSomeconcerns",
                             "b_treat1:overall_robHigh",
                             "sd_study__treat12")) |> 
  mcmc_neff() +
  theme(legend.position = "none" )

Does the Trace-Plot Exhibit Convergence?

Code
mcmc_trace(rob_mod, pars = c("b_treat1",
                             "b_treat1:overall_robSomeconcerns",
                             "b_treat1:overall_robHigh",
                             "sd_study__treat12"))

Does the Histogram Have Enough Information?

Code
mcmc_hist(rob_mod, pars = c("b_treat1",
                             "b_treat1:overall_robSomeconcerns",
                             "b_treat1:overall_robHigh",
                             "sd_study__treat12"))
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Do the Chains Exhibit a Strong Degree of Autocorrelation?

Code
mcmc_acf(rob_mod, pars = c("b_treat1",
                             "b_treat1:overall_robSomeconcerns",
                             "b_treat1:overall_robHigh",
                             "sd_study__treat12"),
         lags = 10)

Does the Posterior Distribution Make Substantive Sense?

Code
mcmc_dens(rob_mod, pars = c("b_treat1",
                             "b_treat1:overall_robSomeconcerns",
                             "b_treat1:overall_robHigh",
                             "sd_study__treat12"),
          transformations = list("b_treat1" = "exp"))

Published vs. Unpublished

Data

Code
filtered_data_publish = 
  data |> 
  filter(N_cef != "NA") |>  # exclude studies with missing data
  rename("death_com" = deaths_com,
         "study" = Study.ID) |> 
  mutate(publish = ifelse(Comparator == "Control", "Unpublished", "Published"))

long_data_publish = 
  filtered_data_publish |> 
  pivot_longer(
   cols = c(death_cef, death_com, N_cef, N_com),
   names_to = c(".value", "group"),
   names_sep = "_"
) |> 
  # important for model fitting
  mutate(study = factor(study),
         treat = ifelse(group == "com", 0, 1) |> factor(levels = c("0", "1")),
         treat12 = ifelse(group == "com", -0.5, 0.5),
         
         publish = factor(publish,
                             levels = c("Published",
                                        "Unpublished"))
         )

Overall Model

Model formula

Code
# Model 4 in DOI: 10.1002/sim.7588
mf = 
  bf(death | trials(N) ~ 1 + study + treat + (treat12 - 1 | study))

Priors

Code
get_prior(mf, family = binomial, data = long_data_publish)
                prior     class                       coef group resp dpar
               (flat)         b                                           
               (flat)         b              studyAl411056                
               (flat)         b              studyAl411070                
               (flat)         b              studyAl411071                
               (flat)         b              studyAl411075                
               (flat)         b              studyAl411079                
               (flat)         b              studyAl411097                
               (flat)         b              studyAl411113                
               (flat)         b              studyAl411118                
               (flat)         b              studyAl411119                
               (flat)         b              studyAl411120                
               (flat)         b              studyAl411137                
               (flat)         b              studyAl411149                
               (flat)         b              studyAl411154                
               (flat)         b              studyAl411157                
               (flat)         b              studyAl411169                
               (flat)         b              studyAl411179                
               (flat)         b              studyAl411184                
               (flat)         b              studyAl411196                
               (flat)         b              studyAl411205                
               (flat)         b              studyAl411206                
               (flat)         b              studyAl411213                
               (flat)         b              studyAl411219                
               (flat)         b              studyAl411221                
               (flat)         b              studyAl411222                
               (flat)         b              studyAl411228                
               (flat)         b              studyAl411230                
               (flat)         b              studyAl411242                
               (flat)         b              studyAl411245                
               (flat)         b              studyAl411247                
               (flat)         b              studyAoun1997                
               (flat)         b              studyAtay2004                
               (flat)         b           studyAufiero1997                
               (flat)         b     studyAvilesMRobles2019                
               (flat)         b           studyBarckow1993                
               (flat)         b         studyBeaucaire1999                
               (flat)         b             studyBiron1998                
               (flat)         b             studyBohme1998                
               (flat)         b          studyBonfitto1999                
               (flat)         b               studyBow2006                
               (flat)         b           studyBradley2019                
               (flat)         b             studyBreen1997                
               (flat)         b         studyCannavino2015                
               (flat)         b      studyChandrasekar2000                
               (flat)         b             studyChang1998                
               (flat)         b            studyCherif2004                
               (flat)         b            studyChuang2002                
               (flat)         b           studyCordero2001                
               (flat)         b        studyCordonnier1997                
               (flat)         b           studyCornely2001                
               (flat)         b           studyCornely2002                
               (flat)         b            studyCPM0896003                
               (flat)         b            studyCPM4497002                
               (flat)         b            studyCPM6796007                
               (flat)         b         studyEdelstein1991                
               (flat)         b             studyErman2001                
               (flat)         b            studyFujita2016                
               (flat)         b            studyGentry1991                
               (flat)         b            studyGentry1992                
               (flat)         b           studyGhalaut2007                
               (flat)         b           studyGlauser1997                
               (flat)         b             studyGomez2001                
               (flat)         b             studyGomez2010                
               (flat)         b          studyGrossman1999                
               (flat)         b         studyHoepelman1993                
               (flat)         b          studyHolloway1996                
               (flat)         b             studyHuang2002                
               (flat)         b             studyHuang2005                
               (flat)         b              studyJehn1998                
               (flat)         b             studyJiang2003                
               (flat)         b            studyJindal2016                
               (flat)         b            studyKebudi2001                
               (flat)         b             studyKieft1994                
               (flat)         b            studyKutluk2004                
               (flat)         b              studyKwon2008                
               (flat)         b               studyLee2018                
               (flat)         b         studyLeophonte1993                
               (flat)         b               studyLin2001                
               (flat)         b               studyLiu2019                
               (flat)         b            studyMallet1997                
               (flat)         b           studyMcCabe1996a                
               (flat)         b           studyMcCabe1996b                
               (flat)         b           studyMustafa2001                
               (flat)         b           studyNakane2015a                
               (flat)         b           studyNakane2015b                
               (flat)         b            studyNaseem2011                
               (flat)         b            studyNewton1993                
               (flat)         b              studyOguz2006                
               (flat)         b                studyOi2020                
               (flat)         b             studyORyan1996                
               (flat)         b          studyPaladino2007                
               (flat)         b     studyPonceMdeMLeon1999                
               (flat)         b           studyPreheim1995                
               (flat)         b              studyQian2023                
               (flat)         b              studyRaad2003                
               (flat)         b           studyRamphal1997                
               (flat)         b      studySaezMLlorens1995                
               (flat)         b      studySaezMLlorens2001                
               (flat)         b            studySaito1992a                
               (flat)         b            studySaito1992b                
               (flat)         b              studySanz2002                
               (flat)         b         studySarashina2014                
               (flat)         b            studySchaad1998                
               (flat)         b           studySchrank1995                
               (flat)         b          studySchwartz1996                
               (flat)         b              studySeo2017a                
               (flat)         b              studySeo2017b                
               (flat)         b studyShahid2008DCPM4495001                
               (flat)         b           studySharifi1996                
               (flat)         b            studyTalaie2008                
               (flat)         b            studyTamura2002                
               (flat)         b             studyUygun2009                
               (flat)         b              studyWang1999                
               (flat)         b            studyWillis1998                
               (flat)         b          studyYakovlev2006                
               (flat)         b           studyZanetti2003                
               (flat)         b            studyZervos1998                
               (flat)         b                     treat1                
 student_t(3, 0, 2.5) Intercept                                           
 student_t(3, 0, 2.5)        sd                                           
 student_t(3, 0, 2.5)        sd                            study          
 student_t(3, 0, 2.5)        sd                    treat12 study          
 nlpar lb ub       source
                  default
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
                  default
        0         default
        0    (vectorized)
        0    (vectorized)
Code
informative = bayesmeta::TurnerEtAlPrior("all-cause mortality",
                                         "pharma",
                                         "pharma")


logmean = informative$parameters["tau", "mu"]
logsd = informative$parameters["tau", "sigma"]

# logmean
# -2.09

# logsd
# 0.705

priors_overall = 
  # study's baseline (treat = 0) log odds
   prior(normal(0, 1.5), class = Intercept) + 
  prior(normal(0, 1.5), class = b) + 
  # overall treatment effect, log odds ratio, "theta"
  prior(normal(0, 0.82), class = b, coef = "treat1") + 
  # between-study SD, "tau"
  prior(lognormal(-2.09, 0.705), class = sd)

Fit the model only sampling from prior distributions

Code
publish_overall_mod_prior = 
  brm(
    sample_prior = "only",
    formula = mf,
    family = binomial,
    prior = priors_overall, 
    data = long_data_publish,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "publish_overall_mod_prior.Rds"),
    file_refit = "on_change"
)

Now check if the estimand distribution looks reasonable (prior predictive check):

Code
as_draws_df(publish_overall_mod_prior) |> 
  ggplot() +
  aes(x = b_treat1) +
  stat_halfeye(.width = 0.95,
               point_interval = median_hdi) + 
  scale_x_continuous(
    breaks = c(0.1, 0.2, 0.5, 1, 2, 5, 10) |> log(),
    labels = function(x) exp(x)
    ) +
  coord_cartesian(x = c(0.1, 10) |> log())  +
  labs(x = "\nOdds Ratio (log scale)",
       y = " ",
       title = "Prior Predictive Check of the Treatment Effect Estimand") +
  theme(axis.text.y = element_blank())

Fit the full model (data + prior)

Code
publish_overall_mod = 
  brm(
    formula = mf,
    family = binomial,
    prior = priors_overall, 
    data = long_data_publish,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "publish_overall_mod.Rds"),
    file_refit = "on_change"
)

Model Diagnostics

Rhat

Code
rhat(publish_overall_mod$fit) |> 
  mcmc_rhat_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios of effective sample size to total sample size

Values are colored using different shades (lighter is better). The chosen thresholds are somewhat arbitrary, but can be useful guidelines in practice.

light: between 0.5 and 1 (high)

mid: between 0.1 and 0.5 (good)

dark: below 0.1 (low)

Code
neff_ratio(publish_overall_mod) |> 
  mcmc_neff_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Does the Trace-Plot Exhibit Convergence?

Code
mcmc_trace(publish_overall_mod, pars = c("b_treat1","sd_study__treat12"))

Does the Histogram Have Enough Information?

Code
mcmc_hist(publish_overall_mod, pars = c("b_treat1","sd_study__treat12"))
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Do the Chains Exhibit a Strong Degree of Autocorrelation?

Code
mcmc_acf(publish_overall_mod, pars = c("b_treat1","sd_study__treat12"),
         lags = 10)

Does the Posterior Distribution Make Substantive Sense?

Code
mcmc_dens(publish_overall_mod,
          pars = c("b_treat1","sd_study__treat12"),
          transformations = list("b_treat1" = "exp"))

Meta-Regressions

Model formula

Code
# Formula with meta-regression
mf_publish = 
  bf(death | trials(N) ~ 1 + study + treat*publish + (treat12 - 1 | study))

get_prior(mf_publish, data = long_data_publish, family = binomial)
                prior     class                       coef group resp dpar
               (flat)         b                                           
               (flat)         b         publishUnpublished                
               (flat)         b              studyAl411056                
               (flat)         b              studyAl411070                
               (flat)         b              studyAl411071                
               (flat)         b              studyAl411075                
               (flat)         b              studyAl411079                
               (flat)         b              studyAl411097                
               (flat)         b              studyAl411113                
               (flat)         b              studyAl411118                
               (flat)         b              studyAl411119                
               (flat)         b              studyAl411120                
               (flat)         b              studyAl411137                
               (flat)         b              studyAl411149                
               (flat)         b              studyAl411154                
               (flat)         b              studyAl411157                
               (flat)         b              studyAl411169                
               (flat)         b              studyAl411179                
               (flat)         b              studyAl411184                
               (flat)         b              studyAl411196                
               (flat)         b              studyAl411205                
               (flat)         b              studyAl411206                
               (flat)         b              studyAl411213                
               (flat)         b              studyAl411219                
               (flat)         b              studyAl411221                
               (flat)         b              studyAl411222                
               (flat)         b              studyAl411228                
               (flat)         b              studyAl411230                
               (flat)         b              studyAl411242                
               (flat)         b              studyAl411245                
               (flat)         b              studyAl411247                
               (flat)         b              studyAoun1997                
               (flat)         b              studyAtay2004                
               (flat)         b           studyAufiero1997                
               (flat)         b     studyAvilesMRobles2019                
               (flat)         b           studyBarckow1993                
               (flat)         b         studyBeaucaire1999                
               (flat)         b             studyBiron1998                
               (flat)         b             studyBohme1998                
               (flat)         b          studyBonfitto1999                
               (flat)         b               studyBow2006                
               (flat)         b           studyBradley2019                
               (flat)         b             studyBreen1997                
               (flat)         b         studyCannavino2015                
               (flat)         b      studyChandrasekar2000                
               (flat)         b             studyChang1998                
               (flat)         b            studyCherif2004                
               (flat)         b            studyChuang2002                
               (flat)         b           studyCordero2001                
               (flat)         b        studyCordonnier1997                
               (flat)         b           studyCornely2001                
               (flat)         b           studyCornely2002                
               (flat)         b            studyCPM0896003                
               (flat)         b            studyCPM4497002                
               (flat)         b            studyCPM6796007                
               (flat)         b         studyEdelstein1991                
               (flat)         b             studyErman2001                
               (flat)         b            studyFujita2016                
               (flat)         b            studyGentry1991                
               (flat)         b            studyGentry1992                
               (flat)         b           studyGhalaut2007                
               (flat)         b           studyGlauser1997                
               (flat)         b             studyGomez2001                
               (flat)         b             studyGomez2010                
               (flat)         b          studyGrossman1999                
               (flat)         b         studyHoepelman1993                
               (flat)         b          studyHolloway1996                
               (flat)         b             studyHuang2002                
               (flat)         b             studyHuang2005                
               (flat)         b              studyJehn1998                
               (flat)         b             studyJiang2003                
               (flat)         b            studyJindal2016                
               (flat)         b            studyKebudi2001                
               (flat)         b             studyKieft1994                
               (flat)         b            studyKutluk2004                
               (flat)         b              studyKwon2008                
               (flat)         b               studyLee2018                
               (flat)         b         studyLeophonte1993                
               (flat)         b               studyLin2001                
               (flat)         b               studyLiu2019                
               (flat)         b            studyMallet1997                
               (flat)         b           studyMcCabe1996a                
               (flat)         b           studyMcCabe1996b                
               (flat)         b           studyMustafa2001                
               (flat)         b           studyNakane2015a                
               (flat)         b           studyNakane2015b                
               (flat)         b            studyNaseem2011                
               (flat)         b            studyNewton1993                
               (flat)         b              studyOguz2006                
               (flat)         b                studyOi2020                
               (flat)         b             studyORyan1996                
               (flat)         b          studyPaladino2007                
               (flat)         b     studyPonceMdeMLeon1999                
               (flat)         b           studyPreheim1995                
               (flat)         b              studyQian2023                
               (flat)         b              studyRaad2003                
               (flat)         b           studyRamphal1997                
               (flat)         b      studySaezMLlorens1995                
               (flat)         b      studySaezMLlorens2001                
               (flat)         b            studySaito1992a                
               (flat)         b            studySaito1992b                
               (flat)         b              studySanz2002                
               (flat)         b         studySarashina2014                
               (flat)         b            studySchaad1998                
               (flat)         b           studySchrank1995                
               (flat)         b          studySchwartz1996                
               (flat)         b              studySeo2017a                
               (flat)         b              studySeo2017b                
               (flat)         b studyShahid2008DCPM4495001                
               (flat)         b           studySharifi1996                
               (flat)         b            studyTalaie2008                
               (flat)         b            studyTamura2002                
               (flat)         b             studyUygun2009                
               (flat)         b              studyWang1999                
               (flat)         b            studyWillis1998                
               (flat)         b          studyYakovlev2006                
               (flat)         b           studyZanetti2003                
               (flat)         b            studyZervos1998                
               (flat)         b                     treat1                
               (flat)         b  treat1:publishUnpublished                
 student_t(3, 0, 2.5) Intercept                                           
 student_t(3, 0, 2.5)        sd                                           
 student_t(3, 0, 2.5)        sd                            study          
 student_t(3, 0, 2.5)        sd                    treat12 study          
 nlpar lb ub       source
                  default
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
             (vectorized)
                  default
        0         default
        0    (vectorized)
        0    (vectorized)

Priors

Code
priors_regression_publish = 
  # study's baseline (treat = 0) log odds, when subgroup == "Published"
  prior(normal(0, 1.5), class = Intercept) +
  prior(normal(0, 1.5), class = b) + 
  # treatment effect in reference subgroup "Published", log odds ratio
  prior(normal(0, 0.82), class = b, coef = "treat1") + 
  # change in log-odds when treat = 0 for subgroup "Unpublished"
  prior(normal(0, 1.2), class = b, coef = "publishUnpublished") + 
  # change in the treatment effect for subgroup "Unpublished" compared to "Published"
  prior(normal(0, 0.3), class = b, coef = "treat1:publishUnpublished") + 
  # between-study SD, "tau"
  prior(lognormal(-2.09, 0.705), class = sd)

Fit the model only sampling from prior distributions

Code
publish_subgroup_mod_prior = 
  brm(
    sample_prior = "only",
    formula = mf_publish,
    family = binomial,
    prior = priors_regression_publish, 
    data = long_data_publish,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "publish_subgroup_mod_prior.Rds"),
    file_refit = "on_change"
  )

Now check if estimand distributions look reasonable (prior predictive check):

Code
as_draws_df(publish_subgroup_mod_prior) |> 
  reframe("Published" = b_treat1,
          "Unpublished" = b_treat1 + `b_treat1:publishUnpublished`) |> 
  pivot_longer(everything()) |> 
  mutate(name = factor(name, levels = c("Published",                     
                                        "Unpublished") |> rev()
  )
  ) |> 
  
  ggplot() +
  aes(x = value, y = name, fill = name) +
  stat_halfeye(.width = 0.95,
               point_interval = median_hdi) +
  scale_fill_manual(values=met.brewer("Isfahan1", 2) |> rev()) +
  scale_x_continuous(
    breaks = c(0.1, 0.2, 0.5, 1, 2, 5, 10) |> log(),
    labels = function(x) exp(x)
  ) +
  coord_cartesian(x = c(0.1, 10) |> log())  +
  labs(x = "\nOdds Ratio (log scale)",
       y = " ") +
  theme(legend.position = "none",
        axis.text.y = element_text(hjust = 1))

Fit the full model (data + prior)

Code
publish_subgroup_mod = 
  brm(
    formula = mf_publish,
    family = binomial,
    prior = priors_regression_publish, 
    data = long_data_publish,
    
    control = list(adapt_delta = .95),
    backend = "cmdstanr", # faster
    cores = parallel::detectCores(),
    chains = 4,
    warmup = 5000, 
    iter = 10000, 
    seed = 123,
    refresh = 0, # omit sampling text
    file = here("models", "publish_subgroup_mod.Rds"),
    file_refit = "on_change"
  )

Model Diagnostics

Rhat

Code
rhat(publish_subgroup_mod$fit) |> 
  mcmc_rhat_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios of effective sample size to total sample size

Values are colored using different shades (lighter is better). The chosen thresholds are somewhat arbitrary, but can be useful guidelines in practice.

light: between 0.5 and 1 (high)

mid: between 0.1 and 0.5 (good)

dark: below 0.1 (low)

Code
neff_ratio(publish_subgroup_mod) |> 
  mcmc_neff_hist() +
  theme(legend.position = "none" )
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Does the Trace-Plot Exhibit Convergence?

Code
mcmc_trace(publish_subgroup_mod,
           pars = c("b_treat1",
                    "b_treat1:publishUnpublished",
                    "sd_study__treat12"))

Does the Histogram Have Enough Information?

Code
mcmc_hist(publish_subgroup_mod,
           pars = c("b_treat1",
                    "b_treat1:publishUnpublished",
                    "sd_study__treat12"))
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Do the Chains Exhibit a Strong Degree of Autocorrelation?

Code
mcmc_acf(publish_subgroup_mod,
           pars = c("b_treat1",
                    "b_treat1:publishUnpublished",
                    "sd_study__treat12"),
         lags = 10)

Does the Posterior Distribution Make Substantive Sense?

Code
mcmc_dens(publish_subgroup_mod,
           pars = c("b_treat1",
                    "b_treat1:publishUnpublished",
                    "sd_study__treat12"),
          transformations = list("b_treat1" = "exp"))